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P representation techniques, which have been very successful in quantum optics and in other 
fields, are also useful for general bosonic quantum dynamical many-body calculations such as Bose- 
Einstein condensation. We introduce a representation called the gauge P representation which 
greatly widens the range of tractable problems. Our treatment results in an infinite set of possible 
time-evolution equations, depending on arbitrary gauge functions that can be optimized for a given 
quantum system. In some cases, previous methods can give erroneous results, due to the usual 
assumption of vanishing boundary conditions being invalid for those particular systems. Solutions 
are given to this boundary-term problem for all the cases where it is known to occur: two-photon 
absorption and the single-mode laser. We also provide some brief guidelines on how to apply the 
stochastic gauge method to other systems in general, quantify the freedom of choice in the resulting 
equations, and make a comparison to related recent developments. 
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I. INTRODUCTION 

One of the most difficult problems in theoretical 
physics is also conceptually the simplest. How does 
one calculate the dynamical time evolution or even the 
ground state of an interacting many-body quantum sys- 
tem? In essence, this is a natural part of practically any 
comparison of quantum theory with experiment. The 
difficulty is that the Hilbert space of all but the most 
trivial cases can be enormous. This implies that a finite 
computer is needed to to solve problems that can easily 
become nearly infinite in dimensionality, if treated using 
an orthogonal basis expansion. 

In this paper, we formally introduce and give exam- 
ples of techniques for treating general bosonic many-body 
quantum systems, which we call gauge P representations. 
These are an extension of the phase-space method called 
the positive-P representation [1], and have been recently 
used in the context of interacting Bose gases [2, 3]. The 
advantages of the new technique are the following. 

(1) The elimination of certain types of mathematical 
terms known as boundary-term corrections, which have 
caused problems in the positive-P representation for over 
a decade [4-6]. This is the main focus of the present 
paper. 

(2) Greatly reduced sampling error in computations. 
Gauge P representations have been used recently to re- 
duce the sampling error in Kerr oscillator simulations [2]. 

(3) The extension of allowable problems to 'imaginary- 
time' canonical ensemble calculations. These problems 
will be treated elsewhere. 

Related extensions to the positive-P representation — 
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although restricted to the scalar interacting Bose gas 
problem — have also been introduced recently. Different 
procedures have been introduced by Carusotto, Castin, 
and Dalibard [7, 8], and by PHmak, Olsen, and Col- 
lett [9]. These methods implicitly assume the absence of 
boundary term corrections. This paper unifies and sub- 
stantially generalizes all these recent advances. It also 
shows how the gauge method can be used to solve the 
long-standing problem of boundary-term corrections in 
the positive P representation. Comparisons to the other 
methods are given in an Appendix. 

Owing to the work of Wilson [10], and many others 
[11], we know that large Hilbert space problems can often 
be treated using stochastic or Monte Carlo techniques for 
the ground-state, particle masses, and finite-temperature 
correlations. This is the basis for much work in com- 
putational quantum statistical mechanics, and in QCD 
as well. However, Wilson's and other related methods 
are restricted to static or 'imaginary-time' calculations, 
rather than quantum-dynamical problems. 

Methods Hke these that use orthogonal basis sets have 
not proven useful for quantum dynamics; owing to the 
notorious phase problem that occurs when trying to sum 
over families of paths in real-time Feynman path inte- 
grals. For this reason, the many-body quantum time- 
evolution problem is often regarded as inherently insol- 
uble due to its exponential complexity. In fact, it was 
this very problem that motivated the original proposal of 
Feynman [12] to develop quantum computers. In these 
(usually conceptual) devices, the mathematical problem 
is solved by a physical system consisting of evolving 
'qubits' or two-state physical devices. Fortunately, this 
method of doing calculations is not the only one, since 
no large enough quantum computer exists at present[13]. 

Historically, an alternative route is the use of quasi- 
probability representations of the quantum state, which 
either implicitly or explicitly make use of a non orthog- 
onal basis. The term quasi-probability is used because 
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there can be no exact mapping of all quantum statos to a 
classical phase space with a positive distribution [14] that 
also preserves all the marginal probabilities. These meth- 
ods include the Wigner [15] (W), Glauber-Sudarshan (P) 
[16, 17], and Husimi (Q) [18, 19] representations. The 
classical phase-space representations can be classified ac- 
cording to the operator ordering that stochastic moments 
correspond to: the W is symmetrically ordered, the Q 
is anti-normally ordered, while the P representation is 
normally ordered. Apart from numerous laser physics 
and quantum optics calculations, these methods have also 
been used to some extent in quantum statistical mechan- 
ics: for example, the theory of BEC phase fluctuations 
[20]. 

None of these methods result in a stochastic time evo- 
lution with a positive propagator when there are nonlin- 
earities. To achieve this, a better approach is to use a 
non-classical phase space of higher dimension. A com- 
plex higher-dimensional 'R representation' was proposed 
in Glauber's seminal paper on coherent state expansions 
[16]. The first probabilistic method of this type was 
the positive-P representation [1] (+P), which has proved 
capable of performing stochastic time-domain quantum 
calculations in some many-body quantum systems [21]. 
This uses a basis of coherent states that are not orthog- 
onal, thus allowing freedom of choice in the construction 
of the representation. The positive-P representation of 
a quantum state is therefore the most versatile out of a 
large group of quasi-probability distributions developed 
to aid quantum mechanical calculations. It has been suc- 
cessfully applied to mesoscopic systems such as quantum 
solitons [21-23] and the theory of evaporative cooling 
[24], which correctly reproduces the formation of a BEC 
— as observed in experiment [25-27]. 

Quasi-probability distributions of this type are com- 
putationally superior to direct density matrix methods, 
which are susceptible to computational complexity blow- 
up for large Hilbert spaces. Provided certain boundary 
terms vanish, the usual procedure is to generate a Fokker- 
Planck equation (which will vary depending on the dis- 
tribution chosen) from the master equation, and then to 
convert this to a set of stochastic Langevin equations. 
For some simple cases, it may even be possible to ar- 
rive at appealing results directly from the Fokker-Planck 
equation (FPE). The resulting stochastic equations can 
be thought of just as quantum mechanics written in dif- 
ferent variables. They have two main advantages over 
orthogonal basis-state methods, as follows. 

First, the whole quantum dynamics can be written ex- 
actly in terms of a small number of stochastic equations. 
In a one-mode case, there is just one complex variable 
for P and Q and W, and two complex variables for +P. 
Although a simulation requires us to average over many 
realizations of the stochastic process, this is often more 
practical than solving the infinite set of deterministic 
equations required to solve directly for all the elements of 
a density matrix. Such an infinite set may be truncated, 
but this is only a good approximation for a system with 



few particles, and no more than a few modes. 

Second, for a many-mode problem the Hilbert space 
dimension \s N = for the case of n particles dis- 
tributed over M modes. This gives exponential growth 
as a function of the number of modes. However, the num- 
ber of quasi-probability dynamical equations grows only 
linearly with the number of modes, rather than exponen- 
tially in the case of direct methods. Other stochastic 
methods, known as quantum-trajectory methods, can be 
used to reduce the N'^ dimensionality ofa.n NxN density 
matrix problem to that of the A'^-dimensional underlying 
Hilbert space — but this is clearly insufficient to solve the 
complexity problem inherent in the exponential growth 
of the Hilbert space dimension. 

There are, however, some caveats when using these 
distributions. In particular, the vanishing of bound- 
ary terms is an important fundamental issue with quasi- 
probability distributions, and it is this issue that we focus 
on mostly in this paper. To get an overall picture, con- 
sider that once we have a time-evolution problem there 
are five typical requirements that are encountered in de- 
riving stochastic equations for quasi-probability represen- 
tations of many-body systems. These requirements occur 
in closed (unitary evolution) systems, in open systems (in 
general described by a master equation) , or even using a 
distribution to solve for the canonical ensemble in imagi- 
nary time. As such, these requirements are generic to the 
use of stochastic equations with operator representations: 

(1) Positive distribution. A well-behaved positive dis- 
tributions for all quantum states, including especially the 
chosen initial condition, is essential for a general algo- 
rithm. For example, a number state has a highly singu- 
lar P distribution, and a W distribution that is negative 
in some regions of phase space [28], making either distri- 
bution impossible to interpret probabilistically for these 
states. The R distribution is inherently complex. Such 
problems do not occur for the Q or +P representation 

— these are positive, and well-behaved for all quantum 
states [1]. 

(2) Ultraviolet convergence. While normally-ordered 
representations are well behaved at large momentum, 
non-normally-ordered representations of quantum fields 

— such as the Q or W representations — typically face 
the problem of ultraviolet divergence in the limit of large 
momentum cutoff [24]. This means that almost any ob- 
servable quantity will involve the simulation of a (nearly) 
infinitely noisy classical field, leading to diverging stan- 
dard deviations in two or more space dimensions, even 
for linear systems. This rules out the Q and W distribu- 
tions for quantum field simulations in higher than one- 
dimensional environments. 

(3) Second- order derivatives. Only FPEs with sec- 
ond or infinite-order derivatives can be translated into 
stochastic equations [29]. Normally-ordered methods 
such as the P and +P representations can handle most 
commonly occurring nonlinearities and two-body interac- 
tions, with only second-order derivatives. Non-normally 
ordered representations of quantum fields often lead to 
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third- or higher-order partial derivatives in the Fokker- 
Planck equation with no stochastic equivalent. For ex- 
ample, the Wigner representation gives such problems for 
almost any nonlinear term in the master equation. 

(4) Positive-definite diffusion. A Fokker-Planck equa- 
tion must have positive-definite diffusion, to allow simu- 
lation with stochastic processes [29]. When the master 
equation has nonlinear terms, this does not occur with 
any of the classical representations. However, the +P 
representation is guaranteed to always produce positive- 
definite diffusion [1], provided no higher derivative terms 
occur. 

(5) Vanishing boundary terms. In the derivation of 
the Fokker-Planck equations, it is assumed that certain 
boundary terms arising in partial integration can be ne- 
glected. This is not always the case. Boundary terms 
due to power-law tails can occur when there are mov- 
ing singularities that can escape to infinity in finite time. 
In the +P method, such trajectories may cause system- 
atic errors in stochastic averages [6], especially for non- 
integrable dynamical systems. These problems are expo- 
nentially suppressed when linear damping rates are in- 
creased, but can be large at low damping. 

The +P method is often the representation of choice, 
because it satisfies conditions one to four. Gauge repre- 
sentations (G) combined with stochastic methods to be 
treated in this paper, share these advantages with the +P 
representation. However, they can also satisfy the fifth 
requirement — for an appropriate gauge choice — hence 
allowing all of the mathematical problems in simulating 
time evolution to be treated. For this reason, the present 
paper will focus on solving boundarj'-term issues encoun- 
tered with the +P representation for certain nonlinear 
master equations. The overall picture is summarized in 
Table I, as applied to the two-boson anonlinear absorber 
cases treated here in Sec. IV: 

We emphasize that the particular examples treated 
here have a small particle number and extremely low (or 
zero) linear damping. As such, they are soluble using 
other techniques, which allows us to test the accuracy 
of gauge techniques. Our purpose is to demonstrate the 
success of the stochastic gauge method in simple cases 
where boundary terms arise within the +P representa- 
tion. In this way, we can understand more complex situ- 
ations where no exact result is known. 

We will first derive and describe the stochastic gauge 
method in Sees. 11 and III, and subsequently work 
through two examples: First, solving the boundary-value 
problem for the driven one- and two-photon absorber in 
Sec. IV. Second, in Sec. V we will consider the one-mode 
laser at extremely low power, which exhibits boundary 
term errors when very non-optimal starting conditions 
are used. This example will show that gauge methods 
can also be used to remove errors from this system, but 
some judgment must be employed to avoid choosing a 
pathological initial distribution. In the Appendix, we 
compare the methods derived here with recent related 
extensions of the positive-P representation by Carusotto 



and co-workers [7, 8], Plimak et al. [9], and Deuar and 
Drummond [2]. 

Finally, we point out a sixth requirement of contain- 
ing the growth of sampling error: the averages calcu- 
lated from the stochastic Langcvin equations correspond 
to quantum mechanical expectation values only in the 
limit of infinitely many trajectories. Provided boundary 
terms do not occur, the averages will approach the cor- 
rect values — within an acceptable sampling error — for 
sufficiently many trajectories. If this number should in- 
crease rapidly with time, the simulation will only be of 
use for a limited period [2]. 

The problem of growing sampling error can occur even 
when there are no boundary terms, and may be regarded 
as the ultimate frontier in representation theory, just as 
similar issues dominate the theory of classical chaos. This 
is less of a fundamental issue, since the sampling error 
can always be estimated and controlled by increasing the 
number of trajectories. This is simply a matter of moving 
to a clustered, parallel computational model, or repeat- 
ing the calculation many times. Nevertheless, it is of 
great practical significance. The sampHng error problem 
requires careful gauge optimization, and remains an open 
area for investigation. An intelligent choice of gauge can 
often vastly outweigh a brute force computational ap- 
proach, in terms of sampling error. 

II. GAUGE OPERATOR REPRESENTATIONS 

In gauge representations, the density matrix to be com- 
puted is expanded in terms of a coherent state basis. 
For definiteness, we shall focus on the coherent states of 
the harmonic oscillator, which are useful in expanding 
Bose fields; but other choices are clearly possible. The 
expansion kernel is more general than that used in the 
positive-P representation. In order to define the nota- 
tion, we start by introducing a set of boson annihilation 
and creation operators , a|. The operator rii = a\ai is 
therefore the boson number operator for the ith mode or 
site. Boson commutation relations of [ai,a ] = 5ij hold 
for the annihilation and creation operators. 

A. Coherent states 

If a = (ai, . . . , um) is a complex M-dimensional vec- 
tor with ai = Xi -\- ijji, and a = (ai, . . . ,aM) is an M- 
dimensional vector of annihilation operators, then the 
Bargmann coherent state is defined by 

= exp [a • at] |0) = exp [|a|V2] |a) , (1) 

where \a.) is the usual normalized coherent state which 
is a simultaneous eigenstate of all the annihilation opera- 
tors. The inner product of two Bargmann coherent states 
is 

{13* ||a)=exp[a-/3] . (2) 
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Table I: Comparison of phase-space representations as applied to stochastic treatments of a one- and two- boson nonlinear 
absorber. 



Method 


Form of 


UV 


Order of 


Non-negative 


Stochastic 


Boundary term 


Simulated 




Distribution 


converges 


derivatives 


diffusion 


simulations 


removal 


correctly 


W 


Real 


No 


4 


Sometimes 


No 






Q 


Positive 


No 


4 


Yes 


No 






R 


Complex 


Yes 


2 




No 






P 


Singular 


Yes 


2 


No 


No 






+P 


Positive 


Yes 


2 


Yes 


Yes 


No 


Sometimes 


G 


Positive 


Yes 


2 


Yes 


Yes 


Yes 


Yes 



It is important to notice here that \\a) is an analytic 
function of the complex vector a . The following identi- 
ties therefore follow immediately: 



tti II a) = ai II a) 

sni«) = ^11") • 



(3) 



Since ||a) is an analytic function, the notation d/dai 
is interpreted here as an analytic derivative, which can 
be evaluated in either the real or imaginary directions. 



^ M \ ^ II \ • ^ II \ 



(4) 



Since the coherent states are an over-complete basis 
set, any operator can be expanded in more than one way 
using coherent states. For example, the simplest resolu- 
tion of the identity operator is 



7 = 



1 



I a) {cx\d'^'cx 



(5) 



Thus, introducing a second M-dimensional vector /3, 
we can expand any operator O directly as 



O = 



1 



r.2M 



II 



a) {ct\0\f3*){/3*\(f^ad^^f3 



Here, we have introduced 
1 



0{a,f3) 



2M 



{a\0\n ■ 



(6) 



(7) 



B. P representations 

The possibility of expanding any operator in terms of 
coherent states leads to the idea that such an expansion 
can be used to calculate observable properties of a quan- 
tum density matrix p . Historically, this was first pro- 
posed by Glauber and Sudarshan [16, 17], who suggested 



a diagonal expansion of the form 

J P{a)\a) {a\d^^a . 



(8) 



Unlike the direct expansion given above, this has no off- 
diagonal elements. Surprisingly, expansions of this type 
always exist, as long as the function P{a) is defined 
to allow highly singular generalized functions and non- 
positive distributions [28]. 

As these do not have a stochastic interpretation, the 
positive-P representation was introduced [1], which is de- 
fined as 



p = J P^+\a,f3) 



I a) 



'2M , 



(9) 



for an M-mode system. 

It is always possible to obtain an explicitly positive- 
definite distribution of this type [1], with the definition 



p(+)(a,/3) = 



1 



(47r2) 



2\M 



exp 



a + (3* 





a -13* 


2' 




2 





cx + 13* 



(10) 



This form always exists, as do an infinite class of equiv- 
alent positive distributions. Even simpler ways to con- 
struct the positive-P representation are available in some 
cases. For example, if the Glauber-Sudarshan representa- 
tion exists and is positive, then one can simply construct 



P^+\a,fi) = P{ot)5'^^{a-p*) 



(11) 



The stochastic time evolution of the positive-P distri- 
bution does not generally preserve the above compact 
forms, and may allow less compact positive solutions in- 
stead. However, to obtain a time-evolution equation, it 
is necessary to use partial integration, with the assump- 
tion that boundary terms at infinity can be neglected. 
It is these less compact solutions, occurring during time 
evolution with a nonlinear Fokker-Planck equation, that 
lead to power-law tails in the distribution — and hence 
boundary-term problems caused by the violation of the 
assumption that these terms vanish. 
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C. Gauge representations 

A technique for constructing an even more general pos- 
itive distribution is to introduce a quantum complex am- 
plitude fl, which can be used to absorb the quantum 
phase factor. This leads to the result that any Hermitian 
density matrix p can be expanded in an over-complete 
basis ACa), where a" = {ft, a,0), and 



A("a) = n 



\a) {/3* 



{/3*\\a) 

= n||a)(/3*||exp[-a-/3] . (12) 

We define the gauge representation Gilx) as a real, pos- 
itive function that satisfies the following equation: 



p = / G{~a) K{~a) 



a 



G{'a) A("a)-FH.c. 



a 



(13) 



The last line above follows from the fact that jo is a Her- 
mitian density matrix and G(a) is real. Here, H.c. is 
used as an abbreviation for Hermitian conjugate. The 
use of a complex weight in the above gauge representa- 
tion is similar to related methods introduced recently for 
interacting Bose gases [7, 8], except that we multiply the 
weight by a normahzed (positive-P) projector, in order 
to simplify the resulting algebra. 

As an existence theorem that shows that this represen- 
tation always exists, consider the complex solution 



Po(a,/3) = 



1 



t2M 



{a\p\f3*) {(3* \a) 



(14) 



obtained from Eq. 
simply define 



(7), with a phase 9 = arg(Po) , and 



G{'a) = \Po{a,f3)\5'^{n-exp[ie{a,p)]) 



(15) 



In this type of gauge representation, G(a) is a posi- 
tive distribution over a set of Hermitian density-matrix 
elements A-|- A'''. It is simple to verify that, by construc- 
tion 



Tr 



(16) 



For the case of = 1, this representation reduces to 
the positive-P representation, and the kernel ACa ) is a 
projection operator. Since the positive-P representation 
is a complete representation, it follows that another way 
to construct the gauge P representation is always avail- 
able, if one simply defines 



G("a) =p(+)(a,/3)52(n- 1) 



(17) 



As a simple example, a thermal ensemble with no 
bosons per mode gives a diagonal P distribution that is 
Gaussian, so that 



Gt^('a)(xexp - |a|V"o 5"^^ {a.- p*)5'\^-l) . (18) 



One advantage of the proposed representation is that it 
allows more general expansions than the positive-P dis- 
tribution, and also includes the case of the complex P 
representation — which has proved useful in solving for 
non-equilibrium steady-states in quantum systems. 



D. Operator identities 

The utility of these methods arises when they are used 
to calculate time (or imaginary time — for which the 
positive-P distribution cannot be used) evolution of the 
density matrix. This occurs via a Liouville equation of 
generic form 



d_ 

di 



P = Hp) , 



(19) 



where the Liouville superoperator typically involves pre- 
and post-multiplication of p by annihilation and creation 
operators. As an example, the equation for purely uni- 
tary time evolution under a Hamiltonian H is 



H,p 



(20) 



Effects of the annihilation and creation operators on 
the projectors are obtained using the results for the ac- 
tions of operators on the Bargmann states. 



aACa) 
a+A(^) 
A(a) 



[a„ + (3] A(-^) 



(21) 



For brevity, we use d = {dQ,da,dj3) to symbolize ei- 
ther {df = d/dxi) or —i{df = d/dyi) for each of the 
i = 0, ...,2M complex variables a". This is possible 
since A("a ) is an analytic function of 'a, and an explicit 
choice of the derivative will be made later. 

Using the operator identities given above, the operator 
equations can be transformed to an integro-differential 
equation. 



dp 
dt 



^ = / G(7?) /:^A(7?) 



a 



(22) 



Here the anti-normal ordered notation jCa indicates an 
ordering of all the derivative operators to the right. As an 
example, in the Hamiltonian case, if the original Hamilto- 
nian if (a, a^) is normally-ordered (annihilation operators 
to the right), then 



1 



£a = [HA{a, da+f3)- HAil3, dp + a)] 



(23) 



If no terms higher than second order occur, this proce- 
dure gives a differential operator with the following gen- 
eral expansion: 



6 



3 ■ 



(24) 



where, to simplify notation, the Latin indices will 
from now on be summed over i = 1, . . . ,2Af, since no 
derivatives with respect to are used as yet. V is a term 
not involving derivative operators with respect to any of 
the variables in Ix . The drift term A'"^^ that is normally 
found using the positive-P representation is labeled with 
the superscript (+) to identify it. 

At this stage, the usual procedure in representation 
theory is to integrate by parts, provided boundary terms 
vanish. This gives a normally-ordered differential opera- 
tor acting on the distribution itself, of form 



d_ 
dt 



Gia) 



V ■ 



G(^). (25) 



This type of generalized Fokker-Planck equation can 
be treated formally using techniques developed by Gra- 
ham, involving time-symmetric curved-space path inte- 
grals [30]. For computational purposes, we require special 
choices of the analytic derivatives to obtain a positive- 
definite diffusion, so that the path integrals have equiv- 
alent stochastic equations [29]. We emphasize here that 
the equations resulting are quite different to those ob- 
tained from the direct insertion of a coherent state iden- 
tity into a Fcynman path integral — which results in 
severe convergence problems [31]. The usual positive-P 
representation equations are obtained at this stage — 
provided there is no potential term — and can be trans- 
formed to stochastic equations using the techniques de- 
scribed in the following section. 



the non-unique decomposition of the complex diffusion 
matrix D, which determines the stochastic correlations 
in the final equations. Arbitrary functional parameters 
can therefore be inserted into the final stochastic equa- 
tions in the noise coefficients, which may lead to further 
optimization of the simulation. This is because the de- 
composition of the complex diffusion matrix D = BB^ , 
which is needed to define a stochastic process, does not 
specify the resulting noise matrix B completely. 

It has been recently shown by Plimak, Olsen and Col- 
lett [9] that for the Kerr oscillator using a decomposition 
different from the obvious diagonal one leads to impres- 
sive improvements in the signal-to-noise ratio of the sim- 
ulation (briefly described in Appendix A 2). This some- 
what surprising result leads us to try to quantify the 
amount of freedom of choice available from this source. 

Since D = , it can always be diagonalized by a 
complex orthogonal transformation 



(26) 



where A is the diagonal matrix whose square gives the 
eigenvalues of D. Thus 5^+^ = OA can be considered the 
canonical, or "obvious" choice of decomposition, unique 
apart from the 2M signs of the diagonal terms. However, 
for any orthogonal U, if is a valid decomposition 

of D, then so is the matrix B = B^~^^U. Hence, any 
matrix in the whole orthogonal family B = OXU is a 
valid decomposition. This can be easily quantified using 
a basis 

of the M{2M — 1) independent antisymmetric 2M x 2M 
matrices ct^*^) . One simply introduces 



III. GAUGE FUNCTIONS 



iij) 



(27) 



In gauge representations, the time evolution of the rep- 
resentation is modified from the usual positivc-P repre- 
sentation equations, by the introduction of a number of 
arbitrary and freely defined functions on the phase space. 
This freedom of choice is, of course, not present with an 
orthogonal basis, and is due to the non-orthogonal na- 
ture of a coherent basis set. Although we do not investi- 
gate other cases, it is worth noting that a similar gauge 
freedom is implicitly present whenever a non-orthogonal 
expansion is used — even if it involves different states 
from the choice of coherent states made here (e.g., the 
Fock state wave functions in Refs. [7, 8]). 



A. Diffusion gauges 

We first introduce the diffusion gauges, which were im- 
plicitly present in the original positive-P representation, 

but were only recognized recently as allowing improve- 
ments in the sampling error. These gauges occur via 



As an example, for a one-mode case there is one complex 
gauge function introduced this way, which is g'^ — gu- 
The resulting transformation is 



U = exp (5^(12)^ 

= cos(/) + a(i2)sin(5'^) , 



(28) 



where the antisymmetric matrix proportional to 

a Pauli matrix, 



.(12) 



1 
-1 



(29) 



Hence, if the noise was diagonal in the canonical form, 
the transformed (but equivalent) noise matrix becomes 



B = 



All cos(.g'') All sin(.g'') 
-A22 sm{g'^) A22 cos{g'^) 



(30) 
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Now, the 2M-dimensional (complex) orthogonal ma- 
trix family contains M(2M — 1) free complex parame- 
ters, so there are M{2M — 1) diffusion gauge functions 
gijilx ^t) that one can choose arbitrarily. This represents 
a large class of specific gauges that can be used directly 
in simulations, as opposed to the conditions on noise cor- 
relations usually given elsewhere [9]. 

As pointed out by Graham [30], there is a close sim- 
ilarity between the theory of curved-space metrics, and 
path integrals with a space-varying diffusion matrix. In 
the present context, the space is complex, and we have 
a family of gauges that are generated on taking the ma- 
trix square root of the diffusion matrix. We have not 
yet used this matrix square root, but this decomposition 
will be applied to obtain positive-definite equations via 
the choice of analytic derivatives made in the following 
sections. 

The above holds for square noise matrices , but one 
is also free to add more noise coefficients in the manner 
Bq = [Bs,Q]. Then 



BsB^ = D = D-QQ' 



(31) 



and all the 2MW coeflticients in the 2M x W matrix Q 
are additional arbitrary complex functions. The freedom 
in Eg is the same as before [i.e. M{2M — 1) indepen- 
dent complex gauge functions] , with the proviso that Bg 
is now given by OXU where the square of A gives the 
eigenvalues of the modified matrix D. The matrix 
would be unchanged if QQ^ were set to zero, although 
this choice of Q does not appear to be useful; it just adds 
extra noise. In general it is not clear whether or not any 
advantage can be gained by introducing the additional 
off-square gauge functions contained in Q. 

If B is given a functional form dependent on the phase- 
space variables, it may lead to additional terms in the 
Stratonovich form of the equations, which are considered 
later in this section. In this situation one must be careful 
not to introduce additional boundary-term errors arising 
from an excessively rapid growth of the noise gauges. 

There is a subtlety here which one must take some care 
with. The complex noise matrix B is not the matrix that 
usually appears in the theory of stochastic equations. In- 
stead, this matrix is subsequently transformed into an 
'equivalent' stochastic form, by taking advantage of the 
analyticity of the Bargmann states. This means that the 
effect of the diffusion gauges on the final equations also 
makes use of the non-uniqueness of the coherent basis set 
itself. 



derivative term, and also to introduce a stochastic gauge 
to stabilize the resulting drift equations. This defines an 
infinite class of formally equivalent Fokker-Planck equa- 
tions, in a similar way to related procedures in QED and 
QCD. To demonstrate this, we introduce 2M arbitrary 
complex drift gauge functions g = [giCajt)], to give a 
new differential operator Cga whose form differs from 
the original C^J'^ by terms that vanish identically when 
applied to the kernel A("a), 



t-GA = i-A 



+ 



(32) 

The total differential operator Cga has an anti-normal 
Fokker-Planck form. Extending the drift and diffusion 
matrices to include the extra variable we can write this 

— summing repeated a, b, c indices over a = 0, . . . , 2M 

— as 



Cga = 



^ada + ^Dabdadi 



(33) 



The total complex drift vector is A = [Aq, Ai, 
where 

Ao = nv 

Aj = Af^-gkB^k ■ 



(34) 



The new diffusion matrix D. with elements Dab is not 
diagonal, but it can be factorized. Explicitly, it is now a 
square (2M -|- 1) x (2M -|- 1) complex matrix, given by 



D = 



Bg^Q BB^ 



ng 

B 




Og^ B^ 



BB' 



(35) 



Thus, we now have a new stochastic noise matrix with 
one added dimension. 



B = 



ng 

B 



(36) 



The operator (32) was chosen to give this form for B, 
so that the only change in noise is for the O variable. 



B. Drift gauges 

While the diffusion gauges can control sampling er- 
ror due to the correlations of noise terms, they cannot 
ehminate boundary terms due to singular trajectories in 
the drift equations. The extra variable fl allows the dn 
identity to be used to convert any potential term to a 



C. Positive-definite diffusion 

It is always possible to transform these second- 
derivative terms into a positive semi-definite diffusion 
operator on a real space, which is a necessary require- 
ment for a stochastic equation. When D = B^B^ , divide 
B_ = Bf +iB^ into its real and imaginary parts. A similar 

procedure is followed for A . 
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Recalling that the original kernel was analytic, thus al- 
lowing for more than one choice of derivatives, the choice 
for da can now be made definite by choosing it so that 
the resulting drift and diffusion terms are always real, 

Aada - Aldl + Aldl, (37) 
DMb - Bl,Bl,dldt + Bl,Bt,dldt + {x^y). 

Hence, the gauge differential operator can now be written 
explicitly as 



C-GA 



(38) 



where the indices /i, v cover the {AM + 2)-dimensional 
phase-space of the real and imaginary parts of a , so that 
5 = ("?, 'y), and = d/da^. The diffusion matrix D = 

T 

-B B is now positive semi-definite, since, by construction 



Numerical simulations are usually done in the 
Stratonovich calculus, due to superior convergence prop- 
erties [32], so the equivalent complex Stratonovich equa- 
tion allows us to write efficient algorithms. 



daa = dxa + idua 



- 2 (Bbkdb) Bak 



dt + BakdWk , (44) 



where (B^kdi,) = {B^^d^ + B^f,d^). The derivative terms 
above are the Stratonovich correction in the drift, corre- 
sponding to related terms obtained in curved-space path 
integrals. 

These gauge terms are now utilized to stabilize 
coherent-state paths entering into highly non-classical re- 
gions of phase space. This allows one to benefit from 
the over-completeness of coherent states, in reducing the 
sampling error and eliminating boundary terms. 



B = 



^ 

By 



(39) 



so that the diffusion matrix is the square of a real ma- 
trix — explicitly. 



D = 





By 



(40) 



As Cga is now explicitly real as well as positive-definite 
by construction, it can be applied to the Hermitian con- 
jugate kernel as well, resulting in the final time-evolution 
equation. 



dp 

at 



^ = / G(5) CGAHa) 



a 



(41) 



On integrating by parts, provided boundary terms 
vanish, at least one solution will satisfy the following 
(normally-ordered) positive-definite Fokker-Planck equa- 
tion — with the differential operators on the left, each 
acting on all terms to the right. 



dG_ 
'dt 



CgnG 



-d^d^D^^ 



G. 



(42) 



This implies that we have an equivalent set of Ito 
stochastic differential equations available, with 2M real 
Gaussian noises dWi , which are 



dn = n{Vdt + gkdWk 



da 4 



{A^^^ -gkBjk)dt + BjkdWk 



(43) 



The noises obey {dWidWj) = dijdt, and are uncorrelated 
between time steps. 



D. Moments 

The procedure for calculating observable moments is 
slightly different for the gauge representation than for 
the positive P. Any moment can be written in terms of 
the normally ordered operator products and their 

expectation values are given by 



'^stoch 



quant 



(Q + n*) 



(45) 



sl.och 



which differs from the positive-P situation whenever fl 
differs from unity. 

The average norm (17) is always preserved if there is 
no potential term {V = ), since the resulting equation 
for the weight variable is 



dfl = flgkdWk ■ 



(46) 



The decorrelation property of Ito equations [29] then im- 
plies that 



{dn) = {ngk){dWk) = o . 



E. Gauge properties 



(47) 



We turn briefly here to the question of gauge classifica- 
tion and properties. Just as in QED, the over-complete 
nature of the coherent-state expansion means that many 
equivalent, stable gauges exist. However, they may not 
be equivalent in terms of boundary terms. These are de- 
termined by the tails of the distribution function, which 
depends intimately on the gauge chosen for the time evo- 
lution. It is essential that the distribution tails are suf- 
ficiently bounded to eliminate boundary terms arising in 
partial integration. It is sufficient to bound tails better 
than any inverse power law, for which it is conjectured to 
require (as a necessary condition) that all deterministic 
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trajectories are bounded over any finite time interval [6]. 
This issue is discussed in greater detail below, and in Ref. 
[33]. 

The main criteria for a useful gauge are the elimination 
of boundary terms and the reduction of sampling error. 
However, there is an enlarged space of variables for the 
Fokker-Planck equation here. For this reason, it is possi- 
ble to stabilize trajectories in the usual positive-P phase 
space, while introducing new gauge-induced boundary 
terms in the f2 space. When it comes to the formation 
of boundary terms, the phase of CI is generally innocuous 
provided the gauge is periodic in this variable, but the 
gauge distribution must be strongly bounded as |f2| ^ oo 
to prevent new boundary terms from arising. 

We can classify gauges according to their real or imagi- 
nary nature, and their functional dependence; which can 
be on just the phase-space variables, just the quantum 
phase, or on both. This gives rise to nine gauge types, 
depending on the following criteria. 



1. Gauge complexity 

Gauges are in general complex functions, which leads 
to the following classification of gauge complexity: 

1. Real gauge 

2. Imaginary gauge 

3. Complex gauge 

In general, we find that trajectories can be stabiHzed by 
real, imaginary or complex gauges, provided they have 
some (a, 13) phase-space dependence. 

It is worthwhile to note that the imaginary and real 
parts of the gauges affect the behavior of sampling er- 
ror differently. In the Ito calculus, the evolution of the 
weight O due to the gauges is simply dCt = CtgkdWk- 
Typically, i.e., when there are no significant correlations 
between the phase of a ( or 0) and f2, the weight factor 
appearing in moment calculations is just approximately 
Re[f7]. As a general rule, sampling errors are partially 
due to stochastic fluctuations in the phase-space trajec- 
tories, and partially due to stochastic fluctuations in the 
weight function. Thus there is a trade-off; a gauge that 
is strongly stabiHzing may reduce phase-space fluctua- 
tions at the expense of increased weight variance, and 
vice versa. 

To understand the different types of gauges in some- 
what greater detail, we consider the evolution of the 
weight variance for real and imaginary gauges, in a sim- 
ple case where gauge and weight are decorrelated, with 
Q = 1 initially. Let Cl = fl' + ifl" and Qk = g'k + w'k > 
then 

dn' = (n'g'k-n"g'^)dWk , 

dn" = {n'g',^ + n"g'^)dWk . (48) 



If we consider the evolution of the squares of these terms, 
the Ito rules of stochastic calculus give 

dunf) = m'g',-n"g'if)dt, 

d{[n"f) = {{n'g'l + n"g'^f)dt . (49) 

Suppose for simplicity that the gk and CI are approxi- 
mately uncorrelated, then we have two cases to consider. 

1. Real gauge: 

d{[Cl'f) = ([fif)dT, (50) 

where d,T = {gkgk)dt. This initiafly leads to Hnear 
growth in the variance, and hence in the sampHng 
error. The real part of the gauge will cause noise 
directly in CI' , producing asymetric spreading in CI' , 
which can lead to a few rare very highly weighted 
trajectories for times r > 1. The effect of the real 
gauge may become misleading once the distribu- 
tion becomes highly skewed, as the rare trajectories 
that are important for moment calculations may be 
missed if the sample is too small. At long times, 
if (gkgk) is constant and uncorrelated with Ct, then 
the growth becomes exponential, with ([f^']^) = e"^. 

2. Imaginary gauge: 

d{[Cl']^) = {[Cl"f)dT, 
d{[Cl"]^) = ([Of )(ir . (51) 

where dr = {g'lg'^)dt. This leads initiafly to quad- 
ratic growth in the variance of CI' , and hence a 
slower growth in the sampling error. If {gkgk) is 
constant and remains uncorrelated with CI, then 
the growth is given by ([Of) = cosh(T), ([O'f ) = 
sinh(r). An imaginary gauge will cause mutual 
canceling of trajectories that have weights of ran- 
domly positive and negative sign once r > tt. This 
can also have deleterious effects for smafl samples, 
if the average sample weight becomes negative — of 
course, this cannot be true over the entire stochas- 
tic population. 

The generic behavior is more complex than in the exam- 
ples given above, due to correlations between the gauge 
and the normalization. 

Clearly any type of gauge tends to cause growth in the 
norm variance. However, there is an exception to this 
rule: the norm-preserving gauges. This class of gauges is 
of special interest as they generate trajectories having an 
invariant normalization, so that Re[rfO] = 0. From the 
equation for the norm variance, Eq. (49), it follows that a 
necessary and sufficient condition for a norm-preserving 
gauge is that Cl'g'f. = Cl"g'^ . If O' = 1 initially, this 
implies that gk = iCl* fk = «(1 — iCl")fk , where fk is 
a real function. Unless gk = 0, norm-preserving gauges 
are generally functions of both the phase-space variables 
and the weight O . A preliminary study of these gauges 
has shown that these gauges can greatly reduce sampflng 
error, although gauge-induced boundary terms are also 
possible [2] , depending on the choice of fk ■ 
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2. Functional dependence 

From the above analysis, we see that gauges can func- 
tionally depend on any phase-space variable, as well as 
the generalized quantum phase variable or weight f2 . 
This leads to three functional types: 

1. Autonomous (depends on fi only) 

2. Space dependent (depends on phase-space only) 

3. Mixed (depends on aU components of "a including 

n ) 

Autonomous gauges appear to be the least useful since 
they do not affect a or /3 behavior, but gauges of either 
purely space-dependent or mixed type can be used. 

A possible caveat with mixed gauges is that they may 
be much harder to analyze, as two-way couplings will 
occur between the normal phase-space variables a, /3 and 
the weight. 

IV. NONLINEAR ABSORBER CASE 

The nonlinear absorber is an example of a nonlinear 
master equation that can give either correct or incorrect 
results when treated with the usual positive-P represen- 
tation methods, if the boundary terms are ignored. Gen- 
erally, problems only arise when the linear damping has 
exceptionally small values or the number of bosons per 
mode is small (see Fig. 2), so this is not a practical prob- 
lem in optics. However, for other physical systems such 
as a BEG this may be significant. It is a well-studied 
case, and a detailed treatment can be found in Ref. [6]. 
It also has the merit that exact solutions can be readily 
found using other means. By analyzing this example we 
can ensure that the modifications to the drift equations 
obtained from gauge terms, do eliminate boundary terms 
and give correct results. 

Consider a cavity mode driven by coherent radiation, 
and damped by a zero temperature bath that causes both 
one and two photon losses. We have scaled time so that 
the rate of two-photon loss is unity. Without this nonlin- 
ear process, nothing unusual happens. The scaled one- 
photon loss rate is 7, and e is the scaled (complex) driving 
field amplitude. The master equation is 

— = [sa^ — e*a, + —{2apa^ — dp — pa} a) 

+ ^{2a^pa^^ - a^'^a^p - pa^'^a^) . (52) 

Following the treatment of Sec. II, we arrive at the 
gauge representation Stratonovich stochastic equations 

da = [e - a.{al3 + ig + (7 - 1) /2)]dt + iadW , 
df3 = [e* -P{aP + ig +{-/-!) /2)]dt + i/3dW . 
dn = Sndt + n[gdW + gdW] . (53) 



Here Sndt is the appropriate Stratonovich correction 
term [given by the derivative terms in Eq. (44)], which 
depends on the particular gauges chosen. 

With no gauge {g = g = 0), the positive-P 
Stratonovich equations are recovered, 

da = [s - a{aP + - l}/2)]dt + iadW , 

dp = [s* - l3{al3 + {'y-l}/2)]dt + ipdW . (54) 

We will concentrate on the various simplifications of 
this model, which correspond to existing literature, and 
simpler analysis. 

A. Relevance to many-body problems 

The nonlinearity seen here can occur directly in the 
form of a nonlinear collisional damping term in a many- 
body system, so that it can be referred to gcncrically as 
'two-boson absorption'. This type of damping is common 
both to nonhnear photonic and atomic interactions. 

It is of nearly the same form as for an 'imaginary-time' 
thermal equilibrium state calculation for the usual model 
of an alkali-metal Bose gas or BEG [34]. There, for ex- 
ample, the interaction energy between identical bosons of 
mass m and s-wave scattering length in D-dimensional 
space is given by 

H = [ d^x^t2(x)^2(^) ^ (55^ 

m J 

provided that Os is much smaller than other characteristic 
lengths of the system (which is usuahy the case). The 
master equation for an imaginary-time calculation is 

f = 4{^-M^,^e}^, (56) 

where pe is the thermal canonical ensemble density ma- 
trix, p is the chemical potential, is the number oper- 
ator for the entire system, and r = l/fcsT is an inverse 
temperature. Apart from the fact that it is not trace- 
preserving, this is a nonlinearity very similar to that oc- 
curring in the nonlinear absorber master equation. 

While boundary-term discrepancies only occur with 
this nonhnearity for low occupations per mode (see also 
Fig. 2), for a many-mode system at finite temperature 
one expects a large number of modes to have just such 
a low occupation. Thus, it is important to check that 
boundary terms are indeed eliminated. Note that the 
gauge representation simulation is efficient over a wide 
range of occupation numbers. See, for example. Fig. 3. 
More details of applications to both real and imaginary 
time many-body systems with many modes will be given 
elsewhere. 



B. Two-boson absorber 

In its simplest form, corresponding to 7 = e = 0, only 
two-boson absorption takes place. We expect that for a 
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state IV') = J2n '^'il") even boson number components 
will decay to vacuum, and all odd-numbered components 
will decay to |1), leaving a mixture of vacuum and one- 
boson states at long times. 

The positivc-P representation has been found to give 
erroneous results [4, 35-37] due to the existence of mov- 
ing singularities [6], which cause power-law tails in the 
distribution leading to boundary terms. The moment 
usually concentrated on in this system is the number of 
bosons n = a^a, which corresponds to the statistical av- 
erage of n = a/3 in the positive-P representation. This 
has a convenient closed equation (Stratonovich), 



0.55 



dn = —n{n + ig — l/2)rfr + indW"^ 



(57) 



with dW+ = {dW + dW), t = 2t, and g ^ {g + g)/2. 

Let us examine the behavior of the above equation, 
when g = 0, i.e., in the standard, un-gauged formulation. 
The deterministic part of the evolution has a repellor at 
n = 0, and an attractor at n = i. The noise is finite, 
and of standard deviation y/dt/2 at the attractor. We 
can see that the deterministic part of the evolution has 
a single trajectory of measure zero which can escape to 
infinity along the negative real axis. 




2 3 

r (scaled time units) 

Figure 1: Comparison of two-boson damping simulations. 
Circles: positive P simulation; solid line: circular gauge simu- 
lation; dashed line: exact calculation (truncated number-state 
basis). Simulation parameters: 40 000 trajectories; step size 
= 0.005; initial coherent state. Stratonovich semi-implicit 
method [32]. 



a = -(3 = 



(58) 



where tq = l/a{Qf = -l/n(0). This moving singularity 
is known to cause the power law behavior of the Fokker- 
Planck solution at large \n\, which means that integration 
by parts is not in fact valid — which leads to incorrect 
results. 

Indeed, it can be easily seen that in the steady-state 
limit, all trajectories in a simulation will head toward 
n = \, making limt^oo(?^) = \- Quantum mechanics, 
however, predicts that if we start from a state po, the 
steady state will be 



With the gauge (61), the Stratonovich equations be- 
come 



dn 
dQ. 



-n[ n 



l/2)dT indW+ 



(62) 



Q.{[n+{n- \n\f]dT/2 + i{n - \n\)dW+ } 



Phase-space trajectories have changed now, but since it 
has all come from the same master equation, it still de- 
scribes the same system. Consider the equations for the 
polar decomposition of n = re^"^ , 



lim {fi) 



Ed 

j=0 



2j|po|l + 2j) 



For a coherent state |ao) input, say, this will be 



lim (n) = i ("l - e-2|"«l') . 



(59) 



(60) 



Thus we can expect that the positive P simulation will 

I 1^ 

give correct results only when e'""! ^1. 

To correct the problem we have to change the phase- 
space topology in some way to prevent the occurrence of 
moving singularities. We have found that a good gauge 
for a two-boson absorber nonlinearity in general is 



g = g = i{n — \n\) . 



(61) 



This replaces the —n^ term in Eq. (57), which may be- 
come repulsive from zero, with —n\n\ which is always a 
restoring force, and so never leads to super-exponential 
escape. 



dr = — r(r — l/2)dr, 
# = dW+ . 



(63) 



This is exact, and shows that now we have an attractor 
on the circle |n| = i, and a repellor at n = 0, with free 
phase diffusion in the tangential direction. Once trajec- 
tories reach the attractor, only phase diffusion occurs. 
Some more complicated evolution is occurring in the f2 
variable. In any case, there are now no moving singular- 
ities anywhere in the phase space, and simulations corre- 
spond exactly to quantum mechanics. 

Figure 1 compares results for a truncated number-state 
basis calculation, a positive-P calculation, and a "circu- 
lar" gauge (61) calculation for an initial coherent state 
of ao = l/\/2. Figure 2 compares steady-state values 
for exact, positive-P, and gauge calculations for various 
initial coherent states in a wide range. It is seen that the 
gauge calculation is correct to within the small errors due 
to finite sample size. 
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{ 7i(r = 0) ) (number of bosons) 

Figure 2: Steady state expectation values of boson number 
(h.) obtained by gauge simulations (double triangles) com- 
pared to exact analytic results from Eq. (60) (solid line) and 
positive-P simulations (circles) for a wide range of initial co- 
herent states. Size of uncertainty in gauge results due to finite 
sample size is indicated by vertical extent of 'double-triangle' 
symbol. Steady state was observed to have been reached in all 
simulations by r = 7 or earlier (compare with Fig. 1 and 3), 
lience tliis is tlie time for wliicli tlie simulation data is plotted. 
Simulation parameters: 100 000 trajectories; step size = 0.01. 



C. One- and two-boson absorber 

If we now turn on the one-boson decay as well, but 
still do not have any driving, wc expect that all states 
will decay to the vacuum on two time scales 1 and I/7. 
If 7 » 1, nothing interesting happens, however if 7 < 1, 
we should first see a rapid decay to a mixture of vacuum 
and one-boson states due to the two-boson process, and 
then a slow decay of the one-boson state to the vacuum 
on a time scale of r w 2/7. 

In this case the positive-P equations display different 
behavior depending on whether 7 is above or below the 
threshold 7=1. Below threshold, we have an attractor 
at n = (1 — 7)/2, and a repellor at n = 0, while above 
threshold, the attractor is at n = 0, and the repellor 
at n = —(7 — l)/2. In either case, there is a singular 
trajectory along the negative real axis, which can cause 
boundary-term errors. It turns out that the steady state 
calculated this way is erroneous while 7 < 1, and there 
are transient boundary-term errors while 7 < 2 [4]. The 
false steady state below threshold hes at the location of 
the attractor: (1 — 7)/2. 

Let us try to fix this problem using the same circular 
gauge (61) as before. The equation for r is now 

dr = -r(r - [1 - 7]/2)(iT , (64) 

while the 4) and evolution is unchanged. So, above 
threshold we are left with only an attractor at n = 0, 
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Figure 3: Comparison of simulations for system with both 
single- and double-boson damping. Relative strength 7 = 
0.1; Circles: positive-P simulation, solid line: circular gauge 
simulation; dashed line: exact calculation (truncated number- 
state basis). Gauge simulation parameters: 10^ trajectories; 
step size varies from 0.0001 to w 0.006; initial coherent state 
|10) with (n) = 100 bosons. 



while below threshold we have a repellor at n = sur- 
rounded by an attracting circle at r = (1 — 7)/2. This 
phase space again has no moving singularities. 

The results of simulations for the parameter 7 = 0.1 
are shown in Fig. 3. The gauge simulation tracks the 
exact results. We have chosen 7 <C 1 so that a sys- 
tem with two widely differing time scales is tested. The 
circular gauge avoids the false results of the positive-P 
simulation. Note also that the gauge simulation remains 
efficient for a wide range of occupation numbers — from 
(n) ~ 100 ^ 1, where the positive P is also accurate, to 
{h) Ri 0.1 <C 1 where it is totally incorrect. 



D. Driven two-boson absorber 

The other type of situation to consider is when we have 
a driving field as well as two-boson damping. In these 
considerations we have set the one-boson damping rate 
to zero (7 = 0), since this process never causes any of the 
simulation problems anyway, but leaving it out simplifies 
analysis. Failure of the positive-P representation method 
has been found in this limit as well [5], and is evident in 
Fig. 4. The equation for n is no longer stand-alone in this 
case, and we must simulate all three complex variables 
as in Eq. (53), the f2 equation being the same as in the 
undriven case (62). 

A treatment of the singular trajectory problem with 
the same circular gauge (61) leads again to correct re- 
sults, as seen in Fig. 4. 
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t (scaled time units) 

Figure 4: Driven two-boson absorber with e = 0.05. Circles: 
positive P simulation (1000 trajectories); solid Ime: circular 
gauge simulation (10^ trajectories); dashed line: exact calcu- 
lation (truncated number-state basis). Step size At = 0.025. 
Initial vacuum state. 



V. THE SINGLE-MODE LASER 

Let us now consider the second quantum system for 
which systematic errors have been seen with the positive- 
P representation. We will sec that the problem here is 
somewhat different than in the previous case. The dif- 
ference is that for two-boson damping, boundary-term 
errors occur even when we choose an optimal (i.e., com- 
pact) initial distribution to represent our starting state, 
whereas here systematic errors occur only for unreason- 
ably broad initial distributions. Nevertheless, since nor- 
mally it is assumed that the initial condition can be of 
arbitrary breadth it is instructive to investigate how this 
problem can be tackled with stochastic gauge methods. 

We have found that stochastic gauges can be used to 
increase the allowable breadth to include all reasonable 
starting conditions, but once one tries to increase the ini- 
tial spread too much, it becomes unlikely that any gauge 
will remove systematic errors, without introducing too 
much sampling (i.e. random) error instead. 



A. The laser model 

Ito stochastic differential equations for a simple pho- 
tonic or atomic laser model that can be derived from the 
positive-P distribution are [5, 6] 

da = (G — a(3)adT + \/Qdri , 

dp = {G - aP)PdT + ^/Qdr|* (65) 

in appropriate scaled variables, with the complex Gaus- 
sian noise drj obeying {drjdrf) = 2dT. In terms of physical 



parameters, we have 

a = a/Vj\f 

= , (66) 

where r is the scaled time, and A/" ^ 1 is a scaling pa- 
rameter that equals the number of gain atoms in a simple 
photonic laser model. Both G the gain parameter and 
Q > G/Af, the noise parameter, are real and positive. 

Since this time we are again interested in the (scaled) 
boson number (n) = (a/?) = (n)/jV, its evolution can be 
written as a closed equation 

dn = -2{n - a){h - b)dT + 2^/QhdW , (67) 

where now the real Gaussian noise obeys {dW dW) = 
dr, and the deterministic stationary points in the 
Stratonovich calculus are 

a = 1(g + VG' + 2Q) , 

b = 1(g- + . (68) 

We find that the stationary point at a is an attractor, 
and at b we have a repellor. Defining A = 6 — fi, we get 

dA = 2A(A + v/G2 + 2Q) + noise , (69) 

which shows that we again have a singular trajectory 
escaping to infinity in finite time along the negative real 
axis for n < 6. 



B. Initial conditions 

Let us consider the usual case of vacuum initial condi- 
tions. A vacuum can be represented by 

P^+\aJ) = 6{a)50) , (70) 

but also by Gaussian distributions of any variance CTq, 
around the above. 

Note: the distribution of n is non-Gaussian, but has a 
standard deviation of (7n ~ V^o-q in both the real and 
imaginary directions. 

It has been found by Schack and Schenzle [5] that for 
the single-mode laser model, a positive-P simulation of 
pumping from a vacuum will give correct answers if the 
usual i5-function initial condition (70) is used, but will 
have systematic errors if the initial condition used has a 
sufficiently large variance (see Fig. 5). We emphasize here 
that this is not a real problem in practical cases, as the 
variance required to cause systematic errors is typicahy 
extremely large, once the scahng needed to obtain the 
usual (approximate) laser model is taken into account. 
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This can be understood because if we have a suffi- 
ciently broad initial distribution, the region of phase 
space that includes the singular trajectory will be ex- 
plored by the distribution. Even if initially an 'ti \b\, the 
region n < b may be subsequently explored due to the 
presence of the noise terms. 

Apart from the obvious i5-function initial condition, 
one might want to try the canonical distribution of 
Eq. (10), which is a standard positive-P representation 
construction [1]. It will not cause problems as its vari- 
ance is (Tq = l/M, which for any reahstic case will be 
very small (i.e., an ^ Schack and Schenzle discov- 
ered anomalous results when they chose ctq = 1, due to 
an erroneous procedure of scaling the equations — while 
not scaling the canonical initial condition in a. Neverthe- 
less, since any crp is supposed to represent the same state, 
insight into what can be achieved using gauge methods 
is gained if we analyze the systematic errors for such a 
relatively large ao. 



C. Gauge corrections 

The Fokker-Planck equation corresponding to Eq. (65) 



IS 



dP ( d d ~ d'^ ^ 

^ = \ ^[fi - G]a + ^[n - G]f3 + 2Q--^ \ P . 

(72) 

We now introduce gauges using the same method as in 
Sec. II. This leads to the Ito stochastic equations 

da = a{G — n)dT — \/Q{g + ig)dT + \fQdr] , 
dp = p{G - n)dT - ^/Q{g - ig)dT + ^/Qdrf* , 
dfl = n[{g- ig)dr] + {g + ig)dr]* ]/2 . (73) 

It is convenient to define a transformed gauge function 
g, which is also arbitrary, such that 



9 = 
9 = 



{& + 0)9 



(74) 



Changing to h and 6 = ln(f2) variables we obtain the 
Stratonovich equation 



dh = 2h{G-h-g)dT + QdT + 2y/QhdW , 
dQ = -"^dT + SedT + gJ^dW , (75) 



with Ssdt being the appropriate Stratonovich correction 
[given by the derivative terms in Eq. (44) ] for a particular 
gauge function g. 



D. Correcting for the moving singularities 

Consider the deterministic evolution of the real part, 
n^, oih = hx + ihy, 

dn^ = -2hl + 2Gn^ + Q + 2nl- 2nxRe[5] + 2nylm[g] . 

(76) 

The moving singularity is due to the — 2n^ leading term 
for negative values of fix. We now consider criteria for 
choosing the drift gauges as follows. 

(1) It is desirable to keep the gauge terms to a min- 
imum because whenever they act the weights of tra- 
jectories become more randomized — see Sec. IIIEl. 
Thus, let us restrict ourselves to functions g that are 
only nonzero for fia; < 0. 

(2) This immediately leads to another restriction on g: 
To be able to use the efficient numerical algorithms in 
the Stratonovich calculus, we must be able to calculate 
the correction term Sq, which depends on derivatives of 
g\JnlQ. This immediately suggests that g must always 
be continuous, hence, in particular, lim„^_>o( ^ ) = 0. For 
ease of analysis, let us start with a simple form for the 
gauge, g = c — Xhx + ^fiy- This restriction immediately 
implies c = Ay = 0, hence 



9 = 



-ARe[n] if Re[n] < 
if Re[n] > 



(77) 



and Sq = A(Re[n] + n + \n\)/2. when Re[n] < 0, zero 
otherwise. 

(3) The next necessary condition, to remove moving 
singularities, is that the — 2n^ term is canceled, hence: 



A> 1. 



(78) 



(4) Now, if A = 1 there are no systematic errors, but 
the sampling error very quickly obscures everything be- 
cause still heads to — oo exponentially due to the 2Ghx 
term. This takes it into regions of everincreasing \g\, 
and weights quickly become randomized. For slightly 
larger parameters A, the evolution takes trajectories 
to a point lying far into the negative Ux region where 
the two leading terms balance. Here the trajectories sit, 
and quickly accumulate weight noise. It is clear that for 
an optimum simulation all stationary points of fix in the 
nonzero gauge region must be removed. In this system 
this condition is 



(79) 



An example has been plotted in Fig. 5 where we have 
parameters G = 1, Q = 0.25 (leading to a « 1.1124 and 
h w —0.1124 ). We are considering an initial condition of 
(Tq = 0.1, which is already much larger than the canonical 
variance for physically likely parameters. Typical values 
of n initially will be of order an ~ 0.14 > \h\ here. A 
good choice of gauge has A = 4. The use of this gauge 
clearly restores the correct results. 
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T (scaled time units) 

Figure 5: One-mode laser G = 1, Q = 0.25. Dashed line: 
(correct) positive P simulation with delta-tunction initial con- 
ditions (70) o-Q = and 10® trajectories. Dotted- dashed line: 
erroneous positive-P simulation with Gaussian initial condi- 
tions (71) (Jo = 0.1 initially, and 10^ trajectories. Dotted line: 
positive-P simulation with oq = 1, and 10^ trajectories. Solid 
line: gauge calculation for oq =0.1 with A = 4, which corrects 
the systematic error of the positive P. Only 4000 trajectories, 
so as not to obscure other data. Step size in all cases is 0.005. 



E. Non-optimal initial conditions 

As wc increase the spread of the initial distribution 
beyond an ~ |^|, it becomes increasingly difficult to find 
a gauge that will give reasonable simulations. (For ex- 
ample we have tried a wide variety of what seemed like 
promising gauges for CTq = 0.3, with the above values 
of parameters Q and G, and none have come close to 
success). The problem is that while we can remove sys- 
tematic errors, large random noise appears and obscures 
whatever we are trying to calculate. 

Trajectories that start off at a value of n lying signif- 
icantly beyond h require a lot of modification to their 
subsequent evolution to (1) stop them from escaping to 
— oo and (2) move them out of the gauged region of phase 
space so that they do not accumulate excessive weight 
noise. If there are many of these, the trade-off between 
the gauge size and length of time spent in the gauged re- 
gion does not give much benefit anymore. Nevertheless, 
one may be sure that if this is the case, results will at 
worst be noisy and unusable, rather than being system- 
atically incorrect. 

We stress again that this whole matter of non-optimal 
initial conditions is not a major hurdle to dynamical sim- 
ulations because a compact starting distribution is gen- 
erally found very easily. 



VI. CONCLUSIONS 

The positive-P representation is well suited to complex 
quantum mechanical problems, such as many-body sys- 
tems, but has been known for about a decade to have 
systematic errors in some cases of its use — due to non- 
vanishing boundary terms. The gauge P representation, 
a variant on the usual positive-P representation, can be 
used to eliminate boundary terms and consequently all 
the systematic errors that were encountered previously. 
It can also reduce sampling error in a simulation, and al- 
lows 'imaginary time' calculations of thermal equihbrium 
states. The fact that correct results are immediately ob- 
tained in every case where systematic errors were found 
with the positive-P method, is strong evidence that these 
previous problems were indeed due to boundary terms 
caused by moving singularities in the analytically contin- 
ued deterministic equations. Of course, boundary terms 
can occur for other reasons (for example, if the noise term 
grows too rapidly with radius), so caution is still needed 
in the gauge choice. 

The technique appears to be broadly applicable, and 
only requires the recognition of what instabilities in the 
stochastic equations could lead to problems. It does not 
require detailed knowledge of what the boundary terms 
are, provided instabilities are removed. However, we re- 
mark here that the general specification of necessary and 
sufficient conditions to eliminate boundary terms remains 
an open problem, and clearly requires growth restrictions 
on the gauge terms, both in phase space and quantum- 
amplitude space. Care is also required with the choice 
of the gauge and initial distribution. However, using 
unsuitable gauges or initial conditions may only lead to 
large sampling errors, not systematic errors, provided the 
gauge is chosen to eliminate boundary corrections in the 
first place. Sampling error then allows for a confident 
assessment of the magnitude of inaccuracies in a simula- 
tion, which can be supplemented by numerical analysis 
of the distribution tails. 

The main conclusion we come to is that this method 
does, in the cases studied, provide a complete solution 
to the problem of simulation of a many-body quantum 
system in phase space, under conditions where previous 
direct simulation techniques were not practicable. All 
known technical requirements on the path to obtaining a 
stochastically equivalent description to quantum mechan- 
ics, which is applicable to both large and small particle 
numbers, have been satisfied by this method. For this 
reason, we believe that gauge simulations can be used to 
simulate many quantum systems without systematic er- 
rors when carrying out more difficult calculations, where 
no exact result is known. 

These conclusions must be supplemented by the de- 
tailed study of relevant gauges for particular quantum 
systems. We note, however, that the mathematical tech- 
niques employed here for generating stochastic gauges, 
may well be useful for other representations as well as 
the gauge P representation described here. 
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Appendix A: OTHER EXTENSIONS OF THE 
POSITIVE-P REPRESENTATION 



1. The work of Carusotto, Castin, and Dalibard 

Recently, Carusotto, Castin, and Dalibard [7, 8] 
(CCD) have made related extensions to the positive-P 
representation. These were derived for the particular case 
of an interacting scalar Rose gas, and led to a number of 
conditions for an Ito stochastic evolution to be equivalent 
to a master equation evolution. 

It can be shown quite simply that the equations (43) 
generated by the gauge P representation for this Hamil- 
tonian satisfy the CCD conditions. Wc conjecture that 
these provide the most general possible solution to the 
stochastic problem posed by these authors. In particu- 
lar, db = Il[gkdWk-N{4)idB^ + (f>'^dBi)], using the above 
paper's formalism. Our methods can also treat a much 
larger class of Hamiltonians and master equations than 
considered in the CCD treatment. 

In Ref. [7] systematic errors due to boundary terms 
were not considered. However, evolutions satisfying "ex- 
actness" conditions derived using the same procedure can 
contain such errors. 

As an example, following the CCD procedure [7] for a 
one-mode two-boson absorber master equation, as in Eq. 
(52) with 7 = 5 = 0, one arrives at the conditions 



6^ 



dBidB^ = 
dBf = 

Fi = -dbdBi/U , 

F2 = -db*dB2/n* , 

/ = n(iv</.i<^;)2 , 



(Al) 



(A2) 



where (referring back to the notation in this present pa- 
per). 



#1 = da/y/N = F^dt + dBi 

d(j,2= dp*/VN= F2di + dB2 



dn= d[ne 



-0102^ 



] = fdt + db . 



(A3) 



It can be seen that the positive-P equations (54) satisfy 
these conditions, while producing the erroneous evolu- 
tion seen in Fig. 1. In summary, the methods of the 
CCD paper do not obviate the need to choose gauges 
that eliminate boundary terms. 



2. Noise optimization by Plimak, Olsen, and 
Collett 

In Ref. [9], Plimak, Olsen and Collett have found that 
for some systems (the Kerr oscillator H = woa^d -h 
Ka)'^a? /2, in particular), the most obvious (diagonal) 
choice of noise matrix B may not be the optimal one. 

For example, for the above Hamiltonian, one finds that 
the diffusion matrix (in a, /3) variables is 



D = in 



-a 








= BB^ 



(A4) 



Following the procedure in Eq. (30), an equivalent but 
broader choice of noise matrix B can be any of 



B 



'IK 



iacos{g) iasm{g) 
-f3sm{g) Pcos{g) 



(A5) 



with the usual diagonal decomposition given by 5 = 0. 

However, in Ref. [9] it was found that for a positive-P 
simulation, different decompositions with nonzero con- 
stant g gave the lowest sampling error for coherent state 
initial conditions. In their notation, they introduce 
\/A + 1 = — cos((7), and consider the case of real 
A> 1 (i.e., imaginary g ) only. 



3. Stochastic gauges for the Kerr oscillator 

In Ref. [2] , the sampling error in a Kerr oscillator sim- 
ulation — equivalent to a one-mode EEC model, apart 
from linear terms — was reduced substantially by using 
a representation similar to the gauge P representation 
formally introduced here. The basic differences were the 
following. 

1. Instead of a complex gauge SI, a phase factor e'^ 
with a real 9 variable, was used. 

2. The normalization with respect to the behavior of 
9 was carried out explicitly inside the kernel, rather 
than post-simulation in the moments as in Eq. (45). 

This type of representation is a norm-preserving gauge 
P representation, as discussed earlier. A parametrized 
family of gauges led to stable trajectories (as opposed to 
the large sampling error present with a positive-P simu- 
lation). However, some systematic errors were seen due 
to boundary terms. These boundary terms occurred be- 
cause of the stochastic growth of the gauge term in fl 
space, when 9 approached ±7r/2. With the gauge P 
representation introduced in this paper, a wide range of 
gauges do not lead to any systematic errors [33], provided 
gauge growth is controlled. 

We note here that the norm-preserving gauges have the 
property that, in the present notation, g^ = i[l — iSl"]fk 
. However, while the growth of VL' is stabilized, there 
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is growth in the variance of Q." . This means that the ciently in the weight-function space to avoid finite bound- 
function /fe must behave as a decreasing function of fi" ary terms. The detailed requirements and conditions for 
in order to ensure that the distribution is bounded sufS- this type of gauge wih be treated elsewhere. 
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